Error analysis of free probability approximations to the density of states of disordered systems 
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Theoretical studies of localization, anomalous diffusion and ergodicity breaking require solving 
the electronic structure of disordered systems. We use free probability to approximate the ensemble- 
averaged density of states without exact diagonalization. We present an error analysis that quantifies 
the accuracy using a generalized moment expansion, allowing us to distinguish between different 
approximations. We identify an approximation that is accurate to the eighth moment across all noise 
strengths, and contrast this with the perturbation theory and isotropic entanglement theory. 
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Disordered materials have long been of interest for their 
unique physics such as localization [1 , 2J, anomalous dif- 
fusion J3JIH and ergodicity breaking [5[. Their properties 
have been exploited for applications as diverse as quantum 
dots [6 7|, magnetic nanostructures |8[, disordered met- 
als ||9l ITOl , and bulk heterojunction photovoltaics 1ITTHT31 . 
Despite this, theoretical studies are complicated by the need 
to calculate the electronic structure of the respective systems 
in the presence of random atomic nuclear positions. Con- 
ventional electronic structure theories can only be used in 
conjunction with explicit sampling of thermodynamically 
accessible regions of phase space, which make such calcu- 
lations enormously more expensive than usual single-point 
calculations llTil . 

Alternatively, ensemble-averaged quantities may be com- 
puted or approximated using random matrix theory. In par- 
ticular, techniques from free probability theory allow the 
computation of eigenvalues of sums of certain matrices with- 
out rediagonalizing the matrix sums 1151 . While this has 
been proposed as a tool applicable to general random matri- 
ces Q6 1 and has been used for similar purposes in quantum 
chromodynamics IfPTI , we are not aware of any quantification 
of the accuracy of this approximation in practice. We provide 
herein a general framework for quantitatively estimating the 
error in such situations. We find that this allows us to under- 
stand the relative performances of various approximations, 
and furthermore characterize the degree of accuracy system- 
atically in terms of discrepancies in particular moments of 
the probability distribution functions (PDFs). 

Quantifying the error in approximating a PDF using free 
probability. — We propose to quantify the deviation between 
two PDFs using moment expansions. Such expansions are 
widely used to describe corrections to the central limit the- 
orem and deviations from normality, and are often applied 
in the form of Gram-Charlier and Edgeworth series lfl8l[T9l . 
Similarly, deviations from non-Gaussian reference PDFs can 
be quantified using generalized moment expansions. For 



two PDFs w (£) and w (£) with finite cumulants K\, 
and K\, K.2, . . . , and moments U\, pi, . . . and p.\, p2, ■ ■ ■ respec- 
tively, we can define a formal differential operator which 
transforms w into w and is given by | 
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This operator is parameterized completely by the cumulants 
of both distributions. 

The first k for which the cumulants and jty differ then 
allows us to define a degree to which the approximation 
w ps w is valid. Expanding the exponential and using the 
well-known relationship between cumulants and moments 
allows us to state that if the first k — 1 cumulants agree, but 
the A:th cumulants differ, that this is equivalent to specifying 
that 



w (£) = w (?) 
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At this point we make no claim on the convergence of 
the series defined by the expansion of |[TJ, but use it as a 
justification for calculating the error term defined in ||2j. We 
will examine this claim later. 

The free convolution. — We now take the PDFs to be DOSs 
of random matrices. For a random matrix Z, the DOS is 
defined in terms of the eigenvalues j of the M samples 
Zi, . . . , Z„„ . . . , Z M according to 
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The central idea using free probability to calculate approxi- 
mate DOSs is to split the Hamiltonian H = A + B into two 
matrices A and B whose DOSs, p( > and respectively, 
can be determined easily. The eigenvalues of the sum is 
in general not the sum of the eigenvalues; instead, we ap- 
proximate the exact DOS with the free convolution A EB B, 
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i.e. p( H ^ ~ p( AtSB \ a particular kind of "sum" which can be 
calculated without exact diagonalization of H. The moment 
expansion presented above quantifies the error of this ap- 
proximation in terms of the onset of discrepancies between 

the fcth moment of the exact DOS, ?4 / and that for the free 

approximant fit . By definition, the exact moments are 



n = n 



(AH 



(A + B) 
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where (Z) = E Tr (Z) /N denotes the normalized expected 
trace (NET) of the N x N matrix Z. The kth moment can be 
expanded using the (noncommutative) binomial expansion 
of (A + B)' c ; each resulting term will have the form of a 
joint moment (A ni B mi ■ ■ ■ A" r B nr ) with each exponent n s ,m s 
being a positive integer such that Yl s =\ ( n s + wis ) = k. The 
free convolution fa is defined similarly, except that A and 
B are assumed to be freely independent, and therefore that 
each term must obey, by definition 11211 , relations of the form 



(nU(A"-(^))(B ms -D)) (5a) 
A" s B™ 3 ) + lower order terms, (5b) 
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where the degree k is the sum of exponents n s , m s and the 
second equality is formed by expanding the first line us- 
ing linearity of the NET. For k < 3, this is identical to the 
statement of (classical) independence |21J. 

For arbitrary matrices A and B, we can construct a free 
approximant 



Z = A + Q- 1 BQ, 



(6) 



where Q is a N x N random matrix of Haar measure. For 
real symmetric A and B it is sufficient to consider orthogonal 
matrices Q, which can be generated from the QR decom- 
position of a Gaussian orthogonal matrix |22|. (This can be 
generalized readily to unitary and symplectic matrices for 
complex and quaternionic Hamiltonians respectively.) The 
effect of the similarity transformation Q -1 ■ Q is to apply a 
random rotation to the basis of B with respect to A, and so 
in the N — > oo limit of large matrices, the density of states 
jo( z ) converges to the free convolution AEBB |Tl5ll23l , i.e. that 
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V-]c ' — f° r a ^ k. This provides a numerical sampling 

method for calculating the moments of the free convolution. 

Testing for }t[ A+B ^ ^ p\f^ B ^ then reduces to testing 
whether each centered joint moment of the form in | |5a| is 
statistically nonzero. The cyclic permutation invariance of 
the NET means that the enumeration of all the centered joint 
moments of degree k is equivalent to the combinatorial prob- 
lem of generating all binary necklaces of length k, for which 
efficient algorithms exist l24ll . 

The procedure we have described above allows us to as- 
cribe a degree k to the approximation p( H ^ w p( A ® B *) given 
the splitting H = A + B. For each positive integer n, we 
generate all unique centered joint moments of degree n, and 
test if they are statistically nonzero. The lowest such n for 
which there exists at least one such term is the degree of 
approximation k. This is the main result of our paper. We ex- 
pect that k > 4 in most situations, as the first three moments 



of the exact and free PDFs match under very general condi- 
tions [25J. However, we have found examples, as described 
in the next section, where it is possible to do considerably 
better than degree 4. 

Decomposition of the Anderson Hamiltonian. — As an illus- 
tration of the general method, we focus on Hamiltonians of 
the form 



H 



J h 2 

V 



(7) 



/ h N J 



where / is constant and the diagonal elements hj are identi- 
cally and independently distributed (iid) random variables 
with probability density function (PDF) p^ (?). This is a 
real, symmetric tridiagonal matrix with circulant (periodic) 
boundary conditions on a one-dimensional chain. Unless 
otherwise stated, we assume herein that hj are normally dis- 
tributed with mean and variance o~ 2 . We note that a/ J 
gives us a dimensionless order parameter to quantify the 
strength of disorder. 

So far, we have made no restrictions on the decomposition 
scheme H = A + B other than p^) and p^ being easily 
computable. A natural question to pose is whether certain 
choices of decompositions are intrinsically superior to others. 
For the Anderson Hamiltonian, we consider two reasonable 
partitioning schemes: 
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H = A 2 + B 2 = 
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We refer to these as Scheme I and II respectively. For both 
schemes, each fragment matrix on the right hand side has a 
DOS that is easy to determine. In Scheme I, we have p^ = 
Ph since A\ is diagonal with each nonzero matrix element 
being iid. B\ is simply / multiplied by the adjacency matrix 
of a one-dimensional chain, and therefore has eigenvalues 
A„ = 2/ cos (2nn/N) ||26). Then the DOS of B a is p Bl (?) = 
Hn=i $ (£ — h-n) which converges as N — >■ oo to the arcsine 
distribution with PDF p AS (?) = 1/ (n^AJ 2 - £ 2 ) on the 
interval [— 2|/|,2|/|]. In Scheme II, we have that p^ 2 = 

hi 



PB 



px where px is the DOS of X 



/ 



. Since X 



has eigenvalues e± (?) = h x (?) /2 ± Jh\ (?) /4 + J 1 , their 
distribution can be calculated to be 
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Figure 1: Calculation of the DOS, p(£), of the Hamiltonian H of 
with M = 5000 samples of 2000 x 2000 matrices for (a) low, (b) 
moderate and (c) high noise (a/ J=0.1, 1 and 10 respectively with 
u = l). For each figure we show the results of free convolution 
defined in Scheme I (p( A i fflB i); black solid line), Scheme II (pi^Bz). 
green dashed line) and exact diagonalization (p( H >; red dotted line). 
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Figure 2: Diagrammatic expansion of the term 
(A\B\ A\BiA\Bi A\Bi) in terms of allowed paths dictated by 
the matrix elements of A x and Bj of Scheme I in |8a}. 
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Numerical free convolution. — We now calculate the free 
convolution A EB B numerically by sampling the distribu- 
tions of A and B and diagonalizing the free approximant 
j6j. The exact DOS p( A+B ) and free approximant jo( AfflB ) are 
plotted in Figure [lja)-(c) for both schemes for low, moderate 
and high noise regimes (cr/ } =0.1, 1, 10 respectively). We 
observe that for Scheme I we have excellent agreement be- 
tween jo( H ) and p( A i^ B i) across all values of o~ / }, which is 
evident from visual inspection; in contrast, Scheme II shows 
variable quality of fit. 

We can understand the starkly different behaviors of 
the two partitioning schemes using the procedure outlined 
above to analyze the accuracy of the approximations p( H > ~ 
p{A x m x ) and p (H) _ p (A 2 mB 2 )^ For Scheme I, we observe that 
the approximation ^ is of degree k — 8; the discrepancy 
lies solely in the term (^(AiBi) A ^ (27l . Free probability ex- 
pects this term to vanish, since both A\ and B\ are centered 
(i.e. (A\) = (B\) — 0) and hence must satisfy (5b i with 
tii — m \ — ' ' ' — n i — m A — !• I* 1 contrast, we can calculate 
its true value from the definitions of A\ and B\. By defini- 
tion of the NET (•}, only closed paths contribute to the term. 

\4\ 



Hence, only two types of terms can contribute to < (A\B\ 
these are expressed dia grammatically in Figure [2] The ma- 
trix A\ weights each path by a factor of h, while B\ weights 
each path by /, and in addition forces the path to hop to an 
adjacent site. Consequently, we can write explicitly 



(A 1 B 1 
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=2/ 4 E (hi) 2 E (/z?) + 2/ 4 E (hfY = + 2/V 4 , 

(10) 

where the second equality follows from the independence 
of the hj's. As this is the only source of discrepancy at 



the eighth moment, this explains why the agreement be- 
tween the free and exact PDFs is so good, as the leading 
order correction is in the eighth derivative of p( A i^ B i) with 
coefficient 2ct- 4 / 4 /8! = (crjf /20160. In contrast, we ob- 
serve for Scheme II that the leading order correction is at 
k = 4, where the discrepancy lies in (A^B 2 ). Free prob- 
ability expects this to be equal to (A\Bl) = (A 2 ) (B 2 ) = 

(X 2 ) 2 = (/ 2 + c 2 /2) 2 , whereas the exact value of this 
term is / 2 (/ 2 + c 2 ). Therefore the discrepancy is in the 
fourth derivative of p( AEB ) with coefficient (-c" 4 /4) /4! = 
-a 4 /96. 

Analytic free convolution. — Free probability allows us also 
to calculate the limiting distributions of p( AEBB ) in the macro- 
scopic limit of infinite matrix sizes N — >■ oo and infinite 
samples M — > oo. In this limit, the DOS p( A ® B ) is given as 
a particular type of integral convolution of p^ A > and p( B \ 
We now calculate the free convolution analytically in the 
macroscopic limit for the two partitioning schemes discussed 
above, thus sidestepping the cost of sampling and matrix 
diagonalization altogether. 

The key tool to performing the free convolution analyti- 
cally is the R-transform r (w) = g -1 (w) — if -1 [28 [, where 
g^ 1 is defined implicitly via the Cauchy transform 
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(11) 



For freely independent A and B, the R-transforms lin- 
earize the free convolution, i.e. R( AfflB ) (w) = (w) + 
R( B ) (w), and that the PDF can be recovered from the Plemelj- 
Sokhotsky inversion formula by 
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MSB) {w) = R (AmB) {w) + 
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As an example, we apply this to Scheme I with each iid 
hi following a Wigner semicircle distribution with PDF 
Pw (?) = \/4 - £ 2 /47r on the interval [-2,2]. As described 
earlier (Using semicircular noise instead of Gaussian noise 
simplifies the analytic calculation considerably.) From the 
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Figure 3: DOS, p(g), of the Hamiltonian with M = 5000 samples 
of 2000 x 2000 matrices with (a) low, (B) moderate and (c) high 
semicircular on-site noise (c/ 7=0.1, 1 and 10 respectively with a = 
1), as calculated with exact diagonalization (red dotted line), free 
convolution (black solid line), and perturbation theory with A\ as 
reference (blue dashed line) and Bi as reference (gray dash-dotted 
line). The partitioning scheme is Scheme I of (8a) . 
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DOS p( Al > — pw, we calculate its Cauchy transform (i.e. its 
retarded Green function) 



W 40 Jr z - (£ + le) 



z - Vz 2 - 4 



Next, take the functional inverse 



:(Ai) 



(a;) = w + 



(13) 



(14) 



Subtracting 1/w finally yields the R- transform A A ^ (to) = w. 
Similarly with p( Bl ) = p A g, we have its Cauchy transform 



G^) (z) = lim / r l% ^4 , d£ 



1 



P (Bl) (g) 

iZ-(? + ie)" 5 Vz 2 - 47 2 



and its functional inverse 



Vl+4/W 



(15) 



(16) 



which finally yields the K-transform R( Bl ^ (id) = 
(-1 + v/l+47 2 o; 2 ) /a>. 

To perform the free convolution analytically, we add the 
R-transforms to get (w) = R^ (w) + R (Bl) (id), 

from which we obtain 



JAtBBi) 



(w) = w + 



w 



(17) 



The final steps are to calculate the functional inverse 
g(A t mBi)\ and take its i ma gi nar y par t to obtain p( A i m i). 



Unfortunately, yg( Al EBl ' j cannot be written in a compact 
closed form; nevertheless, the inversion can be calculated nu- 
merically. We present calculations of the DOS as a function 
of noise strength o~/ J in Figure |3j showing again that the free 
convolution is an excellent approximation to the exact DOS. 

Comparison with other approximations. — For comparative 
purposes, we also performed calculations using standard 
second-order matrix perturbation theory [29J for both par- 
titioning schemes. The results are also shown in Figure [3] 
Unsurprisingly, perturbation theory produces results that 



vary strongly with a / J, and that the different series, based 
on whether A is considered a perturbation of B or vice versa, 
have different regimes of applicability. Furthermore it is 
clear even from visual inspection that the second moment of 
the DOS calculated using second-order perturbation theory 
is in general incorrect. In contrast, the free convolution pro- 
duces results with a more uniform level of accuracy across 
the entire range of cr/J, and that we have at least the first 
three moments being correct |25J. 

It is also natural to ask what mean-field theory, another 
standard tool, would predict. Interestingly, the limiting 
behavior of Scheme I as N — > oo is equivalent to a form 
of mean-field theory known as the coherent potential ap- 
proximation (CPA) [30— 32J in condensed matter physics, 
and is equivalent to the Blue's function formalism in quan- 
tum chromodynamics for calculating one-particle irreducible 
self -energies |17|. The breakdown in the CPA in the term 

^ (AiBi) 4 ^ is known ITI I331 ; however, to our knowledge, the 

magnitude of the deviation was not explained. In contrast, 
our error analysis framework affords us such a quantitative 
explanation. 

discuss the predictions of isotropic 
theory, which proposes a linear 
between the classical convolution 

Coo p (A) (0 i° (B) ( x - dx and the free 

convolution p( Am ) (£) in the fourth cumulant |251|34|. The 
classical convolution can be calculated directly from the 
random matrices A and B; by diagonalizing the matrices 
as A — Q^A A Q A and B = Qg 1 A B Q B , the classical convo- 
lution p< A * B ) (£) can be computed from the eigenvalues of 
random matrices of the form Z c ; = A^ + FI AgFI where 
FI is a N x N random permutation matrix. It is instructive 
to compare this with the free convolution, which can be 
sampled from matrices of the form Z' = A^ + Q _1 AgQ, 
which can be shown by orthogonal invariance of the Haar 
measure random matrices Q to be equivalent to sampling 
matrices of the form Z=A + Q _1 BQ described previously. 

As discussed previously, the lowest three moments of Z 
and H are identical; this turns out to be true also for Z c ; [25 J. 
Therefore IE proposes to interpolate via the fourth cumulant, 
with interpolation parameter p defined as 



Finally, we 
entanglement 
interpolation 
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(18) 



We observe that for Scheme I, IE appears to always favor 
the free convolution limit (p = 0) as opposed to the classical 
limit (p = 1); this is not surprising as we know from our pre- 
vious analysis that >c< = %[ lESBl , and that the agreement 
with the exact diagonalization result is excellent regardless 
of c/ J. In Scheme II, however, we observe the unexpected 
result that p is always negative and that the agreement varies 
with the noise strength cr/J. From the moment expansion 
we understand why; we have that the first three moments 



match while k 



(A 2 +B 2 ) (A 2 



-cr 4 /4. The discrepancy 



lies in the term (A^B^), which is expected to have the value 
(A\) (B 2 ) = (7 2 + cr 2 /2) 2 in free probability but instead 
has the exact value J 1 (j 1 + a 1 ). Furthermore, we have that 
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(A 2 



where the only discrepancy lies is in the 



so-called departing term (A2B2A2B2) [25. 34|. This term con- 
tributes to Kf m) but has value (A 2 2 ) (b\) = (j 2 + a 1 /!) 1 
, since for the classical convolution we have that 



m k 



{A 2 *B 2 ) 



(Tl r s=1 (A^B^)) = (Ap^Wsp 1 ™ 5 ). This therefore 



explains why we observe a negative p, as this calculation 
shows that 



^2+52) JA 2 mB 2 ) 
JA 2 *B 2 ) (A 2 m 2 ) 



-2 2 
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(19) 



which is manifestly negative. 

In conclusion, we have demonstrated that the accuracy 
of approximations using the free convolution depend cru- 
cially on the particular choice of partitioning scheme for the 
Hamiltonian. We have found an unexpectedly accurate ap- 
proximation for the DOS of disordered Hamiltonians, both 
for finite dimensional systems and in the macroscopic limit 
N — >■ 00. In particular, this approximation remains accurate 
no matter the strength of noise present in the system. Our 
error analysis framework provides an explanation for this 
accuracy, namely that the lowest seven moments of the eigen- 
values distribution are correct, with the first discrepancy only 
in one particular term arising at the eighth moment. 

We expect our results to be generally applicable to arbi- 
trary Hamiltonians, and are currently investigating the valid- 
ity of these approximations for electronic structure models 
on two- and three-dimensional lattices. These results pave 
the way toward constructing even more accurate approxi- 
mations using free probability, guided by a rigorous error 
analysis framework in terms of the accuracy of successive 
moments. Our results represent an optimistic beginning to 
the use of powerful and highly accurate nonperturbative 
methods for studying the electronic properties of disordered 
condensed matter systems regardless of the strength of noise 
present. We expect that these methods will be especially use- 
ful when the presence of noise is not merely a perturbation 
of a perfect system, but rather, crucial to the emergence of 
unique physical phenomena. 
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